# In this script we use the CIR model and the Vasicek model to estimate the yield
# curve.

# If one wants to clear the R environment
rm(list = ls())

library(tidyverse)

# Load the market yield curve
# (computed on script `3.Estimating the market yield curve by replication.R`)
load("Examples/market.yield.curve.Rdata")

# Check out our data
market.yield.curve %>% rename(`T-t` = t)

########################## CIR model ###################################

# Set parameters

# strength of gravitation back to the average rate
a <- .27

# long-term average interest rate (supposing it is a bit less than the rate at the last period)
b <- .03

# Variance parameter
sigma <- .05

# r(t) at t = 0 (supposing it is a bit more than the rate at the first period)
r_0 <- .052  

(gamma <- sqrt(a^2 + 2 * sigma^2))

# Compute and add to our data all the factors in the affine form of the CIR model
yield.CIR <- market.yield.curve %>%
  rename(`T-t` = t) # ...



    # # Affine term structure
    # A(t,T): ((2*gamma*exp((a+gamma)*`T-t`/2))/((a+gamma)*(exp(gamma*`T-t`)-1)+2*gamma))^(2*a*b/sigma^2),
    # B(t,T): (2*(exp(gamma * `T-t`) - 1))/((gamma + a) * (exp(gamma * `T-t`) - 1 ) + 2 * gamma),
    # P(t,T): `A(t,T)` * exp (-`B(t,T)` * r_0) # Estimated price of €1


 
yield.CIR

# We already have the price of 1€ `P(t,T)`. Let's compute the yield y_t_CIR (in a new column
# called y_t_vsk)
# ...
yield.CIR # ...


# Now let us check how would our immunized portfolio of 3 bonds perform under
# this yield.CIR. Let us load "immuniz.data.Rdata", which was saved in script
# `5.Matching portfolio of bonds and immunization.R`
load("Examples/immuniz.data.Rdata")
immuniz.data

# The PV of the surplus is approx. zero under our market-assumptions-yield
immuniz.data %>% summarize(across(starts_with("PV"),sum))

# Append the three columns of interest from "imuniz.data" to our "yield.CIR" dataset
portfolio.yield.CIR <- yield.CIR #...

####### Portfolio projected behaviour with these approx. yields ######

# Now compute the PV of assets, liabilities and surplus under CIR assumptions
# (PV.CIR.liab, PV.CIR.assets, PV.CIR.surplus)
portfolio.yield.CIR <- portfolio.yield.CIR %>%
  mutate(
    #...
    
    
    )
portfolio.yield.CIR # Not quite zero...

########################### Vasicek model ##############################

# Rescaled sigma parameter (so that the variance term of Vasicek is comparable
# to the one for CIR)
sig.ma <- sigma*sqrt(r_0)

yield.vsk <- market.yield.curve %>%
  rename(`T-t` = t) # ... 

    # # Affine term structure
    # B(t,T) = (1 - exp(-a*`T-t`))/a,
    # A(t,T) = exp(((`B(t,T)` - `T-t`) * (a^2*b - .5*sig.ma^2))/a^2-(sig.ma^2*`B(t,T)`^2)/(4*a)),
    # P(t,T) = `A(t,T)` * exp (-`B(t,T)` * r_0) # Estimated price of €1

yield.vsk

# We already have the price of 1€ `P(t,T)`. Let's compute the yield y_t_CIR
# (in a new column 'y_t_vsk')
# ...

yield.vsk # ...

####### Portfolio projected behaviour with these approx. yields ######

# Now compute the PV of assets, liabilities and surplus under CIR assumptions
# (PV.liab.vsk, PV.assets.vsk, PV.surplus.vsk)
portfolio.yield.CIR <- portfolio.yield.CIR %>%
  mutate(
    #...
    
    
  )
portfolio.yield.CIR # Not quite zero...



###################### Plot the market yield and its approxs. ############

# Construct a data frame with columns `T-t`, y_t.cont, y_t_CIR and y_t_vsk, using
# data frames yield.CIR, yield.vsk.

?left_join
plot.data <- left_join(
  # ...
)

plot.data

# Please put the data in "long format" to facilitate the plot, with a column named
# 'yield type' and 'rate'.
plot.data <- plot.data #...



plot.data

# Now we plot the market yield, the Vasicek approximation and the CIR approximation.

p <- plot.data %>% 
  ggplot(...) # + ... +
  #...
  labs(title = "Market yield and its approximations - Vasicek and CIR") 
p
plotly::ggplotly(p)